#!/usr/bin/perl -w

# Copyright 2014 Thomas H. Schmidt
#
# This file is part of DanceSteps.
#
# DanceSteps is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# any later version.
#
# DanceSteps is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with DanceSteps; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.


### Load Packages ##############################################################
use strict;
use IO::Handle;
use FindBin qw($RealBin); # Absolute path to THIS script.
use lib ($RealBin."/modules/FileIO", $RealBin."/modules");
autoflush STDOUT 1; # For direct output.
 use Math::Trig;

use Commandline;
use PTE;
use FileIO::Basic;
use FileIO::Ifp;
################################################################################



### Default Parameters #########################################################
our $version     = "rc1";             # Version number.
our $year        = "2014";            # Year (of change).

our $verbose     = 0;                 # Be loud and noisy (and global); default: silent.

my $ifpInFile    = 'system.ifp';      # Input GROMOS IFP file.
my $bItpOutFile  = "ffbonded.itp";    # Output GROMACS ffbonded.itp file.
my $nbItpOutFile = "ffnonbonded.itp"; # Output GROMACS ffnonbonded.itp file.

my $helpAndQuit  = 0;                 # Print out program help.
################################################################################



### Internal parameters ########################################################
my %ifpData;                         # Filled by "IFPFiles::readIfp(<IFPFILE>)".
my %itpData;
################################################################################



### Print out program headlines ################################################
printHead();
################################################################################



### Handle commandline parameters ##############################################
addCmdlParam('scalar', 'i',       'Input',       \$ifpInFile,                  $ifpInFile, 'GROMOS file: ifp');
addCmdlParam('scalar', 'ob',      'Output',      \$bItpOutFile,                $bItpOutFile, 'GROMACS file: itp');
addCmdlParam('scalar', 'onb',     'Output',      \$nbItpOutFile,               $nbItpOutFile, 'GROMACS file: itp');
addCmdlParam('flag',   'h',       'bool',        \$helpAndQuit,                $helpAndQuit ? 'yes' : 'no', 'Print help info and quit');
addCmdlParam('flag',   'v',       'bool',        \$verbose,                    $verbose ? 'yes' : 'no', 'Be loud and noisy');

cmdlParser();
################################################################################



### Print program help if the user set the flag ################################
printHelp(getCmdlParamRef(), 1) if $helpAndQuit;
################################################################################



### Read the IFP file ##########################################################
if ($ifpInFile) {
    %ifpData = FileIO::IFP::readIfp($ifpInFile); # Read input IFP file.
    die "ERROR: Cannot find IFP data.\n" unless %ifpData;
}
else {
    printHelp();
}
################################################################################



################################################################################
### Convert bonded parameters ##################################################
################################################################################

### Section: bond-stretching parameters ########################################
for (my $i=0; $i<@{$ifpData{'BONDSTRETCHTYPECODE'}}; $i++) {
    next unless $ifpData{'BONDSTRETCHTYPECODE'}[$i];
    push(@{$itpData{'bonds'}}, sprintf("#define gb_%-2d       %8.6f  %8.4e\n; %s\n;", $i+1, $ifpData{'BONDSTRETCHTYPECODE'}[$i]{'B0'}, $ifpData{'BONDSTRETCHTYPECODE'}[$i]{'CB'}, $ifpData{'BONDSTRETCHTYPECODE'}[$i]{'remarks'}));
}
################################################################################


### Section: bond-angle bending parameters #####################################
for (my $i=0; $i<@{$ifpData{'BONDANGLEBENDTYPECODE'}}; $i++) {
    next unless $ifpData{'BONDANGLEBENDTYPECODE'}[$i];
    push(@{$itpData{'angles'}}, sprintf("#define ga_%-2d       %6.2f      %6.2f\n; %s\n;", $i+1, $ifpData{'BONDANGLEBENDTYPECODE'}[$i]{'T0'}, $ifpData{'BONDANGLEBENDTYPECODE'}[$i]{'CT'}, $ifpData{'BONDANGLEBENDTYPECODE'}[$i]{'remarks'}));
}
################################################################################


### Section: improper (harmonic) dihedral angle parameters #####################
for (my $i=0; $i<@{$ifpData{'IMPDIHEDRALTYPECODE'}}; $i++) {
    next unless $ifpData{'IMPDIHEDRALTYPECODE'}[$i];
    push(@{$itpData{'impdih'}}, sprintf("#define gi_%-2d    %8.5f   %8.5f\n; %s\n;", $i+1, $ifpData{'IMPDIHEDRALTYPECODE'}[$i]{'Q0'}, (rad2deg(sqrt($ifpData{'IMPDIHEDRALTYPECODE'}[$i]{'CQ'})))**2, $ifpData{'IMPDIHEDRALTYPECODE'}[$i]{'remarks'}));
}
################################################################################



### Section: dihedral torsional angle parameters ###############################
for (my $i=0; $i<@{$ifpData{'TORSDIHEDRALTYPECODE'}}; $i++) {
    next unless $ifpData{'TORSDIHEDRALTYPECODE'}[$i];
    push(@{$itpData{'tordih'}}, sprintf("#define gd_%-2d  %8.3f   %8.2f          %d\n; %s\n;", $i+1, $ifpData{'TORSDIHEDRALTYPECODE'}[$i]{'PD'}, $ifpData{'TORSDIHEDRALTYPECODE'}[$i]{'CP'}, $ifpData{'TORSDIHEDRALTYPECODE'}[$i]{'NP'}, $ifpData{'TORSDIHEDRALTYPECODE'}[$i]{'remarks'}));
}
################################################################################


### Write output ITP file ######################################################
backupFile($bItpOutFile) if -e $bItpOutFile;
open(ITPFILE, ">$bItpOutFile") || die "ERROR: Cannot open ITP file \"$bItpOutFile\": $!\n";

print ITPFILE "; Table 2.5.2.1\n;       GROMOS bond-stretching parameters\n;\n;\n;   Bond type code\n;   Force constant\n;   Ideal bond length\n;   Examples of usage in terms of non-bonded atom types\n;\n;\n;   ICB(H)[N]    CB[N] B0[N]\n;\n";
foreach (@{$itpData{'bonds'}}) {
    print ITPFILE $_ . "\n";
}

print ITPFILE ";---\n;       Table 2.5.3.1.\n;       GROMOS bond-angle bending parameters\n;\n;\n; Bond-angle type code\n; Force constant\n; Ideal bond angle\n; Example of usage in terms of non-bonded atom types\n;\n;\n;  ICT(H)[N]  CT[N]  (T0[N])\n;\n";
foreach (@{$itpData{'angles'}}) {
    print ITPFILE $_ . "\n";
}

print ITPFILE ";---\n;       Table 2.5.4.1\n;       GROMOS improper (harmonic) dihedral angle parameters\n;\n;\n; Improper dihedral-angle type code\n; Force constant\n; Ideal improper dihedral angle\n; Example of usage\n;\n;\n; ICQ(H)[N] CQ[N] (Q0[N])\n;\n";
foreach (@{$itpData{'impdih'}}) {
    print ITPFILE $_ . "\n";
}

print ITPFILE ";---\n;       Table 2.5.5.1 (Note: changes with respect to the 43A1 table)\n\n;       GROMOS (trigonometric) dihedral torsional angle parameters\n;\n;\n; Dihedral-angle type code\n; Force constant\n; Phase shift\n; Multiplicity\n; Example of usage in terms of non-bonded atom types\n;\n;\n; ICP(H)[N]  CP[N] PD[N] NP[N]\n;\n";
foreach (@{$itpData{'tordih'}}) {
    print ITPFILE $_ . "\n";
}

close(ITPFILE);
################################################################################

################################################################################
################################################################################
################################################################################



################################################################################
### Convert nonbonded parameters ###############################################
################################################################################

### Section: atomtypes #########################################################
for (my $i=0; $i<@{$ifpData{'SINGLEATOMLJPAIR'}}; $i++) {
#    printf("% 8s %s\n", $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'}, $ifpData{'SINGLEATOMLJPAIR'}[$i]{'matrixLine'});
    my @interactions = split(/\s+/, $ifpData{'SINGLEATOMLJPAIR'}[$i]{'matrixLine'});
    shift(@interactions);
    my $c6Term  = $ifpData{'SINGLEATOMLJPAIR'}[$i]{'sqrtC6'} * $ifpData{'SINGLEATOMLJPAIR'}[$i]{'sqrtC6'};
    my $c12Term = $ifpData{'SINGLEATOMLJPAIR'}[$i]{'sqrtC12_' . $interactions[$i]} * $ifpData{'SINGLEATOMLJPAIR'}[$i]{'sqrtC12_' . $interactions[$i]};

#    print $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'} . "\n";
    my $atomNum = getAtomNum($ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'});
#    print "     => $atomNum\n" if defined $atomNum;

    push(@{$itpData{'atomtypes'}}, sprintf("%5s %4d %10.3f %10.3f %5s %13.10f %13.9e", $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'}, $atomNum, 0, 0, "A", $c6Term, $c12Term));
}
################################################################################



### Section: nonbond_params ####################################################
for (my $i=1; $i<@{$ifpData{'SINGLEATOMLJPAIR'}}; $i++) {
    for (my $j=0; $j<$i; $j++) {
#        print $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'} . "-" . $ifpData{'SINGLEATOMLJPAIR'}[$j]{'atomType'} . "\n";
        my @iInteractions = split(/\s+/, $ifpData{'SINGLEATOMLJPAIR'}[$i]{'matrixLine'});
        my @jInteractions = split(/\s+/, $ifpData{'SINGLEATOMLJPAIR'}[$j]{'matrixLine'});
        shift(@iInteractions);
        shift(@jInteractions);
        my $c6Term  = defined($ifpData{'MIXEDATOMLJPAIR'}[$j][$i]) ? $ifpData{'MIXEDATOMLJPAIR'}[$j][$i]{'c6Term'} : $ifpData{'SINGLEATOMLJPAIR'}[$i]{'sqrtC6'} * $ifpData{'SINGLEATOMLJPAIR'}[$j]{'sqrtC6'};
        my $c12Term = defined($ifpData{'MIXEDATOMLJPAIR'}[$j][$i]) ? $ifpData{'MIXEDATOMLJPAIR'}[$j][$i]{'c12Term'} : $ifpData{'SINGLEATOMLJPAIR'}[$i]{'sqrtC12_' . $iInteractions[$j]} * $ifpData{'SINGLEATOMLJPAIR'}[$j]{'sqrtC12_' . $jInteractions[$i]};

#        push(@{$itpData{'nonbond_params'}}, sprintf("%5s %5s 1   %13.9E   %13.9E", $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'}, $ifpData{'SINGLEATOMLJPAIR'}[$j]{'atomType'}, $c6Term, $c12Term));
        push(@{$itpData{'nonbond_params'}}, sprintf("\t%s\t%s\t1\t%12.6E\t%12.6E", $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'}, $ifpData{'SINGLEATOMLJPAIR'}[$j]{'atomType'}, $c6Term, $c12Term)); # For 1-to-1 comparison with existing 54A7 ff.
    }
}
################################################################################



### Section: pairtypes #########################################################
for (my $i=0; $i<@{$ifpData{'SINGLEATOMLJPAIR'}}; $i++) {
    for (my $j=0; $j<=$i; $j++) {
        my $c6Term  = defined($ifpData{'MIXEDATOMLJPAIR'}[$j][$i]) ? $ifpData{'MIXEDATOMLJPAIR'}[$j][$i]{'c6Term'} : $ifpData{'SINGLEATOMLJPAIR'}[$i]{'lj14pairCS6'} * $ifpData{'SINGLEATOMLJPAIR'}[$j]{'lj14pairCS6'};
        my $c12Term = defined($ifpData{'MIXEDATOMLJPAIR'}[$j][$i]) ? $ifpData{'MIXEDATOMLJPAIR'}[$j][$i]{'c12Term'} : $ifpData{'SINGLEATOMLJPAIR'}[$i]{'lj14pairCS12'} * $ifpData{'SINGLEATOMLJPAIR'}[$j]{'lj14pairCS12'};

#        push(@{$itpData{'pairtypes'}}, sprintf("%5s %5s 1   %13.9E   %13.9E", $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'}, $ifpData{'SINGLEATOMLJPAIR'}[$j]{'atomType'}, $c6Term, $c12Term));
        push(@{$itpData{'pairtypes'}}, sprintf("\t%s\t%s\t1\t%12.6E\t%12.6E", $ifpData{'SINGLEATOMLJPAIR'}[$i]{'atomType'}, $ifpData{'SINGLEATOMLJPAIR'}[$j]{'atomType'}, $c6Term, $c12Term)); # For 1-to-1 comparison with existing 54A7 ff.
    }
}
################################################################################



### Write output ITP file ######################################################
backupFile($nbItpOutFile) if -e $nbItpOutFile;
open(ITPFILE, ">$nbItpOutFile") || die "ERROR: Cannot open ITP file \"$nbItpOutFile\": $!\n";

print ITPFILE "[ atomtypes ]\n;name  at.num    mass     charge ptype            c6             c12\n";
foreach (@{$itpData{'atomtypes'}}) {
    print ITPFILE $_ . "\n";
}

print ITPFILE "\n[ nonbond_params ]\n;   i     j func             c6               c12\n";
foreach (@{$itpData{'nonbond_params'}}) {
    print ITPFILE $_ . "\n";
}

print ITPFILE "\n[ pairtypes ]\n;   i     j func             c6               c12\n";
foreach (@{$itpData{'pairtypes'}}) {
    print ITPFILE $_ . "\n";
}
print ITPFILE "\n";

close(ITPFILE);
################################################################################

################################################################################
################################################################################
################################################################################






################################################################################
### Subroutines ################################################################
################################################################################
sub printHead {
    my @headLines = ("################################################################################",
                     "",
                     "ifp2itp $version",
                     "Convert a GROMOS IFP file to the GROMACS ITP file format.",
                     "Copyright Thomas H. Schmidt, $year",
                     "",
                     "http://code.google.com/p/dancesteps",
                     "",
                     "ifp2itp comes with ABSOLUTELY NO WARRANTY.",
                     "This is free software, and you are welcome to redistribute it",
                     "under certain conditions; type `-copyright' for details.",
                     "",
                     "################################################################################");
    my $maxLength = 80;
    foreach (@headLines) {
        $maxLength = (length $_ > $maxLength) ? length($_) : $maxLength;
    }

    foreach (@headLines) {
        printf "%s%-${maxLength}s\n", ' ' x int(($maxLength - length($_))/2), $_;
    }
}



sub printFoot {
    print <<EndOfFoot;
Please cite:
  [1] Schmidt, T. H. DanceSteps: Dirty toolkit for Molecular Modeling (Manual)
      http://code.google.com/p/dancesteps

EndOfFoot
}



sub printHelp {
    my $cmdLParamRef   = shift;
    my $quitAfterPrint = shift;


    print <<EndOfHelp;
DESCRIPTION
-----------
ifp2itp reads a GROMOS IFP file with van der Waals interactions and converts
it into a GROMACS compatible ITP topology file (ffbonded.itp and ffnonbonded.itp).

USAGE: ifp2itp -i IFPFILE -ob BONDED-ITPFILE -onb NONBONDED-ITPFILE

EndOfHelp

    printParamHelp($cmdLParamRef);

    printFoot();

    exit if $quitAfterPrint;
}



sub printCopyright {
    print <<"EndOfCopyright";
This file is part of DanceSteps.

DanceSteps is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
any later version.

DanceSteps is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with DanceSteps; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

EndOfCopyright
    exit;
}



sub getAtomNum {
    my $atomType = shift;

    if (PTE::isElement($atomType, 1)) {
        return PTE::getElementData($atomType, 'number');
    }
    else {
        $atomType =~ s/[A-Z][a-z]+//g;
#        print " -> " . $atomType . "\n";
        return PTE::getElementData($atomType, 'number') if PTE::isElement($atomType, 1);

        if ($atomType =~ /([A-Z])([A-Z]+)/) {
            my $tempAtomType = $1 . lc($2);
            return PTE::getElementData($tempAtomType, 'number') if PTE::isElement($tempAtomType, 1);
        }

        $atomType =~ s/\d//g;
#        print " -> " . $atomType . "\n";
        return PTE::getElementData($atomType, 'number') if PTE::isElement($atomType, 1);

        if ($atomType =~ /([A-Z])([A-Z])[\-\+]/) {
            $atomType = $1 . lc($2);
        }
#        print " -> " . $atomType . "\n";
        return PTE::getElementData($atomType, 'number') if PTE::isElement($atomType, 1);

        if ($atomType =~ /([A-Z])([A-Z])/) {
            $atomType = $1 . lc($2);
        }
#        print " -> " . $atomType . "\n";
        return PTE::getElementData($atomType, 'number') if PTE::isElement($atomType, 1);

        if ($atomType =~ /([A-Z])[A-Za-z]+$/) {
            $atomType = $1;
        }
#        print " -> " . $atomType . "\n";
        return PTE::getElementData($atomType, 'number') if PTE::isElement($atomType, 1);
    }
    return 0;
}
